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Abstract. A numerical study is presented to analyze the thermal mechanisms of 
unsteady, supersonic granular flow, by means of hydrodynamic simulations of the 
Navier-Stokes granular equations. For this purpose a paradigmatic problem in granular 
dynamics such as the Faraday instability is selected. Two different approaches for the 
Navier-Stokes transport coefficients for granular materials are considered, namely the 
traditional Jenkins-Richman theory for moderately dense quasi-elastic grains, and the 
improved Garzo-Dufty-Lutsko theory for arbitrary inelasticity, which we also present 
here. Both solutions are compared with event-driven simulations of the same system 
under the same conditions, by analyzing the density, the temperature and the velocity 
field. Important differences are found between the two approaches leading to interesting 
implications. In particular, the heat transfer mechanism coupled to the density 
gradient which is a distinctive feature of inelastic granular gases, is responsible for 
a major discrepancy in the temperature field and hence in the diffusion mechanisms. 
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1. Introduction 

The hydrodynamics of granular materials is far from being well understood. The first 
difficulty comes from the kinetic theory level, where the far-from-equilibrium nature of 
the problem leads to both conceptual and technical limitations. Many contributions, 
starting in the '80 of the last century [U [2], have helped to develop a well established 
hydrodynamic theory of granular gases, including mixtures and polydisperse materials. 
However the application to other types of granular materials is still uncertain. 

In academy as well as in industry, one would like to have a good theory for a variety 
of granular flow problems under different conditions. In the process of going from theory 
to real applications, one must resort to good choices of transport coefficients to ensure 
the appropriate modeling of the system. The Navier-Stokes transport coefficients have 
been obtained for dilute and semi-dilute granular gases for selected problems within 
the framework of kinetic theory. However, their validity cannot be guaranteed beyond 
the conditions for which they were derived and as we enter the realm of moderately 
dense materials, where basic assumptions like molecular chaos are not fulfilled. On 
the other hand, a purely empirical approach, like the one used for regular liquids and 
where one measures the transport coefficients, to use them later in the Navier-Stokes 
equations, does not apply for granular hydrodynamics. The reason is that the properties 
of the flow depend strongly and nonlinearly on conditions like the preparation of the 
system, flow rate, and phenomena like dilatancy; plus the fact that, in laboratory 
measurements, effects due to the surface properties of particles, wall roughness, the 
coupling with the interstitial fluid, etc, are generally important. From the theoretical 
point of view, the treatment of granular materials by means of the available statistical- 
mechanics techniques faces inherent difficulties brought out by the dissipative character 
of real grain interactions, which is responsible for microscopic irreversibility, lack of scale 
separation, mesoscopic nature of the flow, and strong nonlinearities in the governing 
equations. 

One of the first attempts to determine the Navier-Stokes transport coefficients 
from the revised Enskog theory was carried out by Jenkins and Richman (JR) [H [2]. 
However, the technical difficulties of the analysis entailed approximations that limited 
their accuracy. In particular, given that their analysis is restricted to nearly elastic 
systems, the inelasticity of collisions only influences the energy balance equation by a 
sink term, and so the expressions of the Navier-Stokes transport coefficients are the same 
as those obtained for elastic collisions. The JR approach has been numerically validated 
by molecular dynamics (MD) simulations in Ref. [3] and in experiments such as granular 
flow past an obstacle [I] and vertically oscillated granular layers [5j El [7J El E]. The 
choice of vibrated granular material as a test case for hydrodynamic theories comes from 
being one of the simplest experiments in which all different regimes of the granular flow 
are present while leading to interesting standing-wave pattern formation and dynamics 
[TU1 ITT] , clustering [12l [13] and phase transitions [TH [151 [16] . 

One of the problems which has been repeatedly revisited as an interesting 
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example of granular collective behavior is the granular Faraday instability. Similarly 
as in that described by Faraday for regular liquids [17], a vibrated granular layer 
develops characteristic patterns of stripes, hexagons and squares, for certain intervals of 
frequencies and amplitudes of oscillation [11] . Experimentally and by means of particle 
simulations, the analogy of vibrated granular materials and regular fluids has been 
clearly shown. In two dimensions, at a certain stage of the motion the Faraday waves 
appear as shown in Fig. [7J while the time evolution of the periodic pattern, with twice 
the periodicity of that of the driving oscillation, can be seen in Fig. |3j 

However, in many applications the dynamics of granular flow is supersonic, in 
a regime where the typical velocity of the flow is often many times or even orders 
of magnitude larger than the thermal velocity. To clarify concepts, the latter is 
proportional to the square root of the fluctuational part of the kinetic energy (which gives 
origin to the so-called granular temperature). More precisely, the vibrational regimes 
which are often imposed in real applications in order to mobilize granular material lead 
to an interplay of alternating diffusion and inertial regimes, giving rise to a rich although 
extremely complex dynamics which we will analyze in detail here. 

Beyond the weak dissipation limit, however, it is expected that the functional 
form of the Navier-Stokes transport coefficients for a granular gas differ from their 
corresponding elastic counterparts. Thus, in subsequent works Garzo and Dufty, and 
Lutsko (GDL) [T8l [19], based on the application of the Chapman-Enskog method 
[20] to the Enskog equation, do not impose any constraints at the level of collisional 
dissipation and take into account the (complete) nonlinear dependence of the Navier- 
Stokes transport coefficients on the coefficient of restitution a. In particular, and in 
contrast to the JR results [U [2] , the heat flux has a contribution proportional to the 
density gradient which defines a new transport coefficient fi, which is not present in the 
elastic case. On the other hand, as for ordinary fluids [20], the Navier-Stokes transport 
coefficients are given in terms of the solutions of a set of coupled linear integral equations 
that are approximately solved by considering the leading terms in a Sonine polynomial 
expansion. In spite of this approximation, the corresponding forms for the transport 
coefficients compare well with computer simulations [2TI [22| 123] , even for quite strong 
inelasticity. 

In a previous paper [7J we studied computationally the Faraday instability 
in vibrated granular disks, comparing the output from particle and Navier-Stokes 
hydrodynamic simulations in detail: the onset of the instability, the characteristic 
wavelength, and the pattern itself by studying the density, temperature and velocity 
fields. This served to validate a Navier-Stokes code for granular material based on 
a WENO (Weighted Essentially Non-Oscillatory) approach [23] which is capable of 
capturing the features of the highly supersonic flow generated by the impact of a piston. 
For this purpose we used the JR expressions [U [2] for the Navier-Stokes transport 
coefficients, valid for elastic hard spheres at moderate densities. The conclusion of the 
study was that the JR results showed qualitative and quantitative agreement with those 
from event-driven MD simulations, in a range of parameters which covered the entire 
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bifurcation diagram of the Faraday instability at the coefficient of restitution a = 0.75. 

As already mentioned, the JR approach however fails describing the heat flux 
accurately, since the transport coefficient /x coupled to the density gradient vanishes in 
the latter approach. The presence of this new term in the heat flux is crucial to explain 
for instance the dependence of the granular temperature with height in MD simulations 
in dilute vibrated systems with gradients only in the vertical direction [25], [261 [27] . Apart 
from that, a value of the coefficient of restitution of 0.75 justifies the use of the correct 
forms of the Navier-Stokes transport coefficients proposed in the GDL approach [181 EH] 
which include the effect of dissipation on momentum and heat transport. 

In the present paper, we follow a similar approach to Ref. [7J, in order to study 
numerically the thermal mechanisms that an oscillating boundary imposes on granular 
material under gravity. That is, we will use the expressions of the Navier-Stokes 
transport coefficients derived in Refs. [TBI HH] to compare the performance of the granular 
Navier-Stokes hydrodynamics with respect to particle simulations. We will also analyze 
the differences between the results provided by the JR approach (TJ [2] and those from 
the current approximation [TBI [T9] to the Navier-Stokes transport coefficients. 

The outline of the paper is as follows: In Section 2 we will review the Navier-Stokes 
theory, introducing the GDL kinetic coefficients for dilute and moderately dense 2D 
granular gases, opposed to the JR kinetic coefficients which are only valid for vanishing 
inelasticity. We will also explain briefly how to treat numerically the Navier-Stokes 
equations, while section 3 will be devoted to the results obtained with JR and GDL 
and their comparison with MD simulations. These will lead to interesting implications 
which will be discussed in more detail in the conclusions section. 

2. Navier-Stokes hydrodynamic theory of granular gases 

We consider a granular fluid composed of smooth inelastic hard disks of mass m and 
diameter a. Collisions are characterized by a (constant) coefficient of normal restitution 
< a < 1. In a kinetic theory description, the relevant information on the system is 
contained in the one-particle velocity distribution function. At moderate densities and 
assuming molecular chaos, the velocity distribution function obeys the (inelastic) Enskog 
kinetic equation [281 129] . Starting from this kinetic theory, one can easily obtain the 
(macroscopic) Navier-Stokes hydrodynamic equations for the number density n(f,t), 
the flow velocity u(f,t), and the local temperature T(f,t) [30] . In the case of two- 
dimensional granular gases, the balance equations read 
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In the above equations, p = mn is the mass density, F is the external force acting on 
the system, P is the pressure tensor, q is the heat flux, and ( is the cooling rate due 
to the energy dissipated in collisions. It is worthwhile to note that the macroscopic 
equations given in Eqs. ([I])-© differ from their counterparts for elastic fluids only via 
the appearance of the cooling rate £ on the right-hand side of Eq. ([3]). On the other 
hand, the corresponding transport coefficients defining the momentum and heat fluxes 
must depend in general on the coefficient of restitution a. 

As it happens for elastic fluids, the usefulness of the balance equations ([I])-® is 
limited unless the fluxes and the cooling rate are specified in terms of the hydrodynamic 
fields and their spatial gradients. To first order in the spatial gradients, the Navier- 
Stokes constitutive equations provide a link between the exact balance equations and 
a closed set of equations for the hydrodynamic fields. The constitutive relation of the 
pressure tensor P^ is 

Pij = pSij - n (djUi + diUj - 5ijV ■ uj - jSyV ■ u, (4) 

where p is the hydrostatic pressure, n is the shear viscosity, and 7 is the bulk viscosity. 
The constitutive equation for the heat flux is 

q = —kVT — /iVn, (5) 

where k is the coefficient of thermal conductivity, and /1 is a new coefficient which does 
not have an analogue for a gas of elastic particles. Finally, to first order in gradients, 
the cooling rate ( can be written as [28] 

C = Co + CiV-m. (6) 

It is important to remark that the derivation of the Navier-Stokes order transport 
coefficients does not limit in principle their application to weak inelasticity. The Navier- 
Stokes hydrodynamic equations themselves may or may not be limited with respect to 
inelasticity, depending on the particular states considered. In particular, the derivation 
of these equations by means of the Chapman-Enskog method assumes that the spatial 
variations of the hydrodynamic fields n, u, and T are small on the scale of the mean free 
path. In the case of ordinary fluids, the strength of the gradients can be controlled by the 
initial or boundary conditions. However, the problem is more complicated for granular 
fluids since in some cases (e.g., steady states such as the simple shear flow jHH [32] ) 
there is an intrinsic relation between dissipation and some hydrodynamic gradient and 
so, the two cannot be chosen independently. Consequently, there are examples for which 
the Navier-Stokes approximation is never valid or is restricted to the quasielastic limit. 
On the other hand, the transport coefficients characterizing the Navier-Stokes order 
hydrodynamic equations are well-defined functions of a, regardless of the applicability 
of those equations. 

As said in the Introduction, the evaluation of the explicit forms of the hydrostatic 
pressure p, the Navier-Stokes transport coefficients 77, 7, k, and fi and the coefficients 
Co and Ci requires to solve the corresponding Enskog equation. However, due to the 
mathematical complexity of this kinetic equation, only approximate results for the above 
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coefficients can be obtained. Here, we consider two independent approaches for hard 
disks proposed by Jenkins and Richman [2] and Garzo and Dufty [18] and Lutsko |19j . 
Let us consider each method separately. 



2.1. Jenkins- Richman (JR) results 

The results derived by Jenkins and Richman [TJ [2] are obtained by solving the Enskog 
equation for spheres [1] and disks [2] by means of Grad's method [33] . The idea behind 
Grad's moment method is to expand the velocity distribution function in a complete 
set of orthogonal polynomials (generalized Hermite polynomials), the coefficients being 
the corresponding velocity moments. Next, the expansion is truncated after a certain 
order k. When this truncated expansion is substituted into the hierarchy of moment 
equations up to order k one gets a closed set of coupled equations. In the case of a 
two-dimensional system, the eight retained moments are the hydrodynamic fields (n, u, 
and T) plus the irreversible momentum and heat fluxes (Pij — pdij and q). 

Although the application of Grad's method to the Enskog equation is not restricted 
to nearly elastic particles, the results derived by Jenkins and Richman [2] neglect the 
cooling effects on temperature due to the cooling rate in the expressions of the transport 
coefficients [see for instance, Eqs. (70), (89), (98), (99), and (100) of Ref. [2] when the 
disks are smooth]. Given that this assumption can only be considered as acceptable for 
nearly elastic systems, the authors of Ref. [2] conclude that their theory only holds in 
the quasielastic limit (a — > 1). 

The explicit forms of the hydrostatic pressure, the Navier-Stokes transport 
coefficients and the cooling rate in the JR theory are given by 

PJR = ^0T[l + (l + a)G(0)], (7) 
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In the above equations, <fi = nix a 1 : /4 is the (dimensionless) volume fraction occupied 
by the granular disks, also called packing fraction, G{4>) = 0x(0), and is ^ ne P & i r 
correlation function. 
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Because of the assumption of near elastic particles in the JR theory, Eqs. (l7|)- (fTT]) 
show clearly that the coefficient of restitution a only enters in the equation of state 
(JTj) and in the expression ( TTTT) for the zeroth-order cooling rate Co- At this level of 
approximation, the expressions of the Navier-Stokes transport coefficients t/jr, 7jr, 
and kjr are the same as those given by the Enskog equation for elastic disks |34j . 

In order to get the dependence of the transport coefficients and the cooling rate 
in both JR and GDL approaches, one has to chose an approximate form for the pair 
correlation function in this paper, we have chosen the forms proposed by Torquato 



x(<f>) 



for < 6 < 



(12) 



for 



< 6 < 



(1-^)2 c -0 

which go through the freezing point <pf = 0.69 and approach the random close packing 
fraction, C = 0.82 with reasonable accuracy. 



2.2. Garzo-Dufty-Lutsko (GDL) results 

The dependence of the Navier-Stokes transport coefficients on the coefficient of 
restitution was first obtained by Garzo and Dufty [18] for hard spheres (d = 3) by 
solving the Enskog equation from the Chapman-Enskog method [20]. These results 
were then extended to an arbitrary number of dimensions by Lutsko [19]. Here, we 
refer to the above theories as the GDL theory. The Chapman-Enskog method [20J is a 
procedure to construct an approximate perturbative solution to the Enskog equation in 
powers of the spatial gradients. As said in the Introduction, the GDL theory considers 
situations where the spatial gradients are sufficiently small and independent of the 
coefficient of restitution a. As a consequence, the corresponding forms of the Navier- 
Stokes transport coefficients are not limited a priori to weak inelasticity since they 
incorporate the complete nonlinear dependence on a. This is the main difference with 
respect to the JR approach. 

On the other hand, as for elastic collisions [20], the Navier-Stokes transport 
coefficients in the Chapman-Enskog method cannot be exactly determined since they are 
defined in terms of the solutions of a set of coupled linear integral equations. It is useful 
to represent these solutions as an expansion in a complete set of polynomials (Sonine 
polynomials) and generate approximations by truncating the expansion. In practice the 
leading terms in these expansions provides an accurate description over the full range 
of dissipation and density since in general they yield good agreement with Monte Carlo 
simulations, except for the heat flux transport coefficients at high dissipation [21] |22) . 
Motivated by this disagreement, a modified version of the first Sonine approximation 
has been recently proposed [23l [36] . The modified Sonine approximation replaces the 
Gaussian weight function (used in the standard Sonine method) by the homogeneous 
cooling state distribution. This new method significantly improves the a-dependence of 
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k and ix since partially eliminates the discrepancies between simulation and theory for 
quite strong dissipation (see for instance, Figs. 1-3 of Ref. [23]). 

The results obtained in the GDL approach for the equation of state and the Navier- 
Stokes transport coefficients for hard disks {d = 2) are 

PGBL = Pm = + (i + 



7T(7 Z 
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where the (reduced) kinetic contributions n* k and ii* k are 
_^ _ 1 + c + |G(0)(1 + a) 2 [2a - 1 + f (1 



a 



(13) 
(14) 
(15) 

(16) 
(17) 

(18) 



2( K - 2Q) 

QK* k (l + (f)d^\nx) + f + |G(0)(1 + a)(l + ^fylnx) [a(a - 1) + ^(14 - 3a + 3a 2 )] 



2< ~ 3Co 

In Eqs. ( JT5"|) - (fT9"l) we have introduced the quantities [36] 
Co = ^x(0)(l-« 2 ) (l + | 

^ = ^x(0)(7-3a)(l + a) (l + g 
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where 



c(a) 



32(l-a)(l-2a 2 



(20) 



(21) 



(22) 



(23) 



57- 25a + 30a 2 (l - a) 
is the fourth cumulant coefficient measuring the deviation of the homogeneous reference 
state from its Gaussian form. Also taking into account Eq. ([12]) . we obtain the expression 
to be used in Eq. (1T91 . 

It is quite apparent that, except the equation of state f|T3|) . the expressions for the 
Navier-Stokes transport coefficients of the GDL results clearly differ from those obtained 
in the JR approach. In fact, Eqs. ( 1T41) . ( 1T51) . ( |T6l) . and ( fTTj) of the GDL theory reduce to 
Eqs. ([8]), ([9]), and (|T0|) . respectively, in the elastic limit (a = 1, and so Co — c = 0). Note 
that the expressions derived by Lutsko [19] neglect in the expressions ( I2T1) and ( 1221) of 
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Figure 1: Bulk viscosity ratio 7GDL/7JR (t°P left), shear viscosity ratio VGBl/VJR 
(top right), thermal conductivity ratio ^gDl/^JR (bottom left), and ^GDL/^JR 
ratio (bottom right) as a function of the restitution coefficient a for three different 
values of the packing fraction <fi: <fi = (solid line), <fi = 0.2 (dashed line), and <fi = 0.4 

(dotted line). 



v* and v*, respectively, the factors of c coming from the non-Gaussian corrections to 
the reference state. These extra factors will be accounted for in our numerical results 
since their effect on transport becomes non negligible at small values of a. In Fig. [T] 
we show the ratio between the bulk viscosity, shear viscosity, and thermal conductivity 
given by the GDL and JR approaches as a function of the coefficient of restitution a 
for different packing fractions <p. Note that the bulk viscosity ratio does not depend on 
<p. We also observe the order of magnitude of the new term in the heat flux due to the 
density gradient in the GDL theory with respect to the heat flux of the JR theory. The 
quantitative percentage of deviation of the transport coefficients with the GDL theory 
from the JR theory is quite significant for a = 0.8 and the different packing fractions </> 
used. We emphasize how the GDL-term related to the density gradient in the heat flux 
becomes very important for a < 0.8. 

Finally, the contributions to the cooling rate are given by 

C ,GDL = £(1 " « 2 ) {^G^) (l + |) , (24) 
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Figure 2: First order correction of the cooling coefficient for GDL theory as a function 
of the coefficient of restitution a for three different values of the packing fraction 0: 
= (solid line), = 0.2 (dashed line), and = 0.4 (dotted line). 



where 



Cl,GDL= oG(0)(l-« 2 
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Equation 021]) agrees with its corresponding counterpart in the JR theory, Eq. (ITT]) , 
when one neglects the non-Gaussian corrections to the reference state (c = 0). Note 
that d vanishes in limits of elastic gases (a = 1, arbitrary volume fraction 0) and of 
dilute inelastic gases (0 = 0, arbitrary values of the coefficient of restitution a). In Fig. 
El we plot the a-dependence of Ci GDL- We observe that the first-order contribution to 
the total cooling rate appears to be more significant as the gas becomes denser. 



2. 3. Numerical scheme for the hydrodynamic granular equations 

The compressible Navier-Stokes-like equations for granular materials ([TJ), (]5]), and 
are solved in conservation form for the convective terms, that is, we numerically solve 
the system for the density, the momentum and the total energy: (n, n«, W) where the 
total energy density W is given by 

1 



W = nT 



-n\u\ 



{21 



This system can be rewritten as a system of nonlinear conservation laws with sources as 
in Ref. [7]. Local eigenvalues and both local left- and right-eigenvectors of the Jacobian 
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matrices of the fluxes are explicitly computable (see Appendix of Ref. [7]). We only 
mention here that the characteristic speeds of the waves in the hyperbolic part of the 
equation can be written in terms of the speed of sound, given by 



2 dp t p dp 



s dn n 2 de ' 

for a general equation of state where p = p(n, e) with the enthalpy e = T for a two 
dimensional system. We refer to Ref. [7] for the full details of the numerical scheme that 
here is applied to both the GDL and the JR Navier-Stokes hydrodynamic equations since 
they share the same structure. Let us just briefly mention that Navier-Stokes terms are 
treated by simple centered high-order explicit in time finite difference approximations 
and considered as sources for the method of lines in the time approximation. Meanwhile 
the Euler (convective) terms are solved in local coordinates by a fifth-order explicit in 
time finite difference characteristic-wise WENO method in a uniform grid following 
Refs. [371 El] • Thus, typical wave speeds and vectors, eigenvalues and eigenvectors of 
the purely hyperbolic part, are correctly resolved. 



3. Results 



We have applied the traditional MD approach to compare the results obtained from the 
Navier-Stokes hydrodynamic equations with the GDL kinetic coefficients, showing also 
those results provided by the previously used JR model as a reference. In all simulations, 
the frequency of the piston motion is f — 3.75 Hz and the amplitude is A = 5.6 particle 
diameters. The system size is tuned to fit three pattern wavelengths in the (horizontal) 
^-direction (125 a), which is periodic. In the (vertical) y-direction, the hydrodynamic 
simulations are constrained into a box of finite height of 60 diameters, whereas the 
MD system is not limited (particles reach the height of 60 diameters very rarely). The 
particles are 783 disks of diameter a = 1 cm and mass m = 1 mg, and g = 9.81 m/s 2 is 
the acceleration of gravity. The coefficient of restitution is a = 0.80, however a similar 
behavior is found regardless the value of the coefficient of restitution between a = 0.60 
and 0.80. At a = 0.85 and beyond instead, the pattern does not form in our system. 
Since JR and GDL will differ less and less, the differences will shrink at high values of a 
anyway. The interesting region is found at intermediate values of a, whereas the use of 
JR is clearly wrong at very low values of the coefficient of restitution. Therefore we will 
show the results for a = 0.80 as a representative case of what one will observe under 
the conditions of the Faraday instability. 

The top and bottom walls in both hydrodynamic simulations are adiabatic and 
impenetrable. More precisely, the normal velocity is zero at the walls, the energy flux 
is zero, and the tangential velocity remains unchanged. The simulation is carried over 
in the comoving frame of the wall, and thus the force per unit mass of the simulated 
system is F = —g(l + Asin(2irft))j, with j = (0, 1). 

We refer the reader to Ref. [7] regarding the details of the averaging procedure 
applied to the MD sequence, here consisting of 1,000 cycles, which leads to the averaged 
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Figure 3: Density field obtained by phase- and space- averaging particle positions from 
the MD simulation. One single wavelength is shown along eight equidistant phases (a) 

to (h). 



(2D) MD hydrodynamic fields for the density (packing fraction) (Fig. |3]), linear 
momentum and thermal energy. From the latter, the temperature field is also obtained. 
These are compared to the corresponding ones generated by the two hydrodynamic 
simulations. 

We disregard the transient originating from the initial condition until the pattern 
of the Faraday instability has fully developed and no changes are observed from period 
to period. After this transient time, which takes about 50 periods of forcing, the system 
reveals a subharmonic periodic dynamics where the period is twice the period of the 
forcing In this regime, we fix the reference time, t — and consider the evolution 
of the profiles of packing fraction, Fig. HJ scaled thermal energy, Fig. [5] scaled granular 
temperature, Fig. |6]and scaled kinetic energy, Fig. [10] using Eq. ( )28|) . The subfigures 
(a-h) correspond to the times t = 0; l /^f~ 1 ] 2 /4/ -1 ; ■ ■ ■] 7 / 4 / _1 - The corresponding 
position of the piston is y = —A sin 2-irft. Despite that the hydrodynamic fields are 
two dimensional, we show ID profiles for a more detailed quantitative analysis. Thus, 
the profiles shown in Figs. HJ|6] are taken at a representative location along the abscissa, 
where the amplitude of the Faraday pattern is developed. The evolution of these profiles 
over the period of excitation is also presented as supplementary online material [38], 
showing the profiles at many more intermediate times. 

3.1. Density 

First of all we are going to discuss the behavior of the packing fraction, Fig. H] Since 
the packing fraction is proportional to the number density <fi = ira 2 n/4, then we will use 
both terms indistinctly. As in subsequent figures, the abscissa represents the height, in 
diameters. On the ordinate we show here the packing fraction. The evolution is shown 
from left to right, and then from top to bottom. Note that the integral of each curve is 
not the same for the hydrodynamics and the MD simulations since it corresponds just 
to a vertical cut at a position where the maximum height of the pattern is achieved. 
Total conservation of mass is maintained in all simulations with high accuracy, see [7] 
for more details. 

At time t = 0, Fig. H(a), the piston is going down through the equilibrium position. 
The height of the material at this location has already grown to a maximum, formed at 



13 




Figure 4: (color online) The profiles of the packing fraction (0) as a function of height 
(in units of er) at selected times over two oscillation periods. For time evolution of the 

profiles see [38] . 
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Figure 5: (color online) Scaled internal energy (<j)T/(ag)) as a function of height (in 
units of a) at selected times over two oscillation periods. For time evolution of the 

profiles see [38] . 
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Figure 6: (color online) Profiles of the temperature (T/(ag)) as a function of height (in 
units of a) at selected times over two oscillation periods for the MD system and the JR 
and GDL solutions. For time evolution of the profiles see |38j . 
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the end of the previous cycle (g, h). Shortly after this time the granular layer experiences 
the impact against the bottom wall and the propagation of a shock wave. Between (a) 
and (c), we see the dissolution of the peak. We observe that the GDL prediction is 
denser than the JR at a distance of 10 diameters from the plate. Just instants following 
frame (c), the layer becomes flat -so does after frame (g), and the material floods to 
neighboring positions to create peaks where valleys previously existed. Shortly after 
(d), another impact with the plate takes place. From frame (d) to frame (g), we see the 
evolution of the density at a valley. 

The MD sequence reveals that the maximum density 0.69 in packing fraction is 
smaller than in both hydrodynamic simulations, reaching the value 0.78. This can 
be due to the irregularity of the MD pattern due to the elasticity of the system 
at a = 0.80, which makes the location of any of the peaks of the MD sequence 
somewhat uncertain. We recall that the granular Navier-Stokes solver does not contain 
fiuctuational -mesoscopic- contributions, while the local noise is enhanced by increasing 
the coefficient of restitution. That is why one needs a factor of 20 times more cycles 
to obtain smooth fields, as compared with the results at a = 0.75, obtained in our 
previous study [7j. There the regularity was much more pronounced, and a much better 
agreement was achieved. 

The role of fluctuating hydrodynamics in granular gases has been an object of 
study for the last decade, since Van Noije and coworkers [3H], or more recently, Brey 
PU] and Costantini |JT]. The essentially mesoscopic dynamics of the granular gas flow 
can not be fully captured by means of macroscopic transport equations. This is easily 
emphasized, for instance, by the need to apply mesoscoping averaging to MD results, in 
order to compare particle and hydrodynamic simulations. Another related effect is the 
diffusion found at the level of the bifurcation threshold of the instability, as observed 
from MD simulations; the hydrodynamic simulations show instead a sharp inception of 
the instability at about T = 2.0, when one represents the wavelength of the Faraday 
pattern as a function of the reduced acceleration T of the plate [7] . 

We expect (minor) differences between the two models in the pattern wavelength, 
however we did not perform the complete analysis of the bifurcation diagram in the 
case of the GDL model for the following reason. As a result of what we have explained 
above, at the threshold of the Faraday instability, the uncertainty in the wavelength as 
determined by the Fourier analysis of the density pattern is quite large [7] . The presence 
of noise turns the transition into a continuous phenomenon, which the hydrodynamic 
simulations without a source of fluctuations cannot exactly reproduce. As a consequence, 
both the JR and the GDL models will be providing somewhat different thresholds, none 
of which will be accounting for the true effect. Beyond the instability threshold, we 
do observe different wavelengths for the case analyzed: 5.6 diameters of amplitude of 
vibration, and reduced acceleration of 2.75. For these parameters, a system 500 wide 
shows 12 wavelengths in the JR simulation, whereas the GDL shows 13. The MD also 
shows 13. This implies that the GDL approach models better this feature as compared 
with JR, at least for the values of the parameters chosen. 
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Figure 7: A snapshot of the MD simulation at the maximum opening of the gap 
(t ~ 0.12 f" 1 ), showing the material stuck at the bottom, between the peaks and the 
plate, whereas there is a completely empty space below the valleys. 



While the GDL and JR profiles do not differ greatly, there are some differences: the 
GDL density is higher at the core of peaks and valleys, as compared to the JR prediction 
at equivalent times. Correspondingly, the packing fraction at the bottom plate is smaller 
in the GDL simulation, and so is the minimum density (0.054 vs. the value of 0.112 
obtained in the JR simulation). However, the minimum density in the averaged MD 
profile is still smaller: 0.004. Also, the impact with the plate occurs later as compared 
with both hydrodynamic simulations, the delay being about 0.16/ -1 . Therefore we may 
argue that in general the accurate expressions of the GDL approximation for the Navier- 
Stokes transport coefficients does not greatly improve the density profile obtained with 
the elastic forms of the JR approach to match the MD results in this problem. A direct 
comparison of the time evolution of densities and velocity fields in full spacial resolution 
can be found in the supplementary material [38] . 

A zoom of the region of the MD system close to the plate during the airborne 
phase will show a few particles stuck to the base of the peaks and empty areas with 
no particles at all below the valleys (Fig. |7J). As a consequence, the impact of the wall 
against the material happens at t — 0.16/ -1 (instead of t = 0). We want to remark that 
this piece of the system is not in the hydrodynamic regime at this moment, but in the 
Knudsen regime, and there is little hope that any hydrodynamic model can reproduce 
this feature in full detail. However the GDL approach to the Navier-Stokes equations 
improves the dynamics of the gap formed as compared with the JR approach in the 
sense that the minimum density at the bottom plate is reduced. On the other hand, 
the density gradients are higher in the GDL theory, a feature which is not observed in 
the MD profiles, which are smoother. The differences are basically due to the presence 
of the coefficient /Uqql (Eq. ( TT9l) of the GDL model), which is absent (/ijR = 0) in the 
elastic case (the JR approach). 
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3.2. Temperature and internal energy 

In Figure [5] we plot the scaled internal energy, <f>T/{ag), where g denotes the gravity 
acceleration. Here we see the evolution of the shock wave travelling across the granular 
layer. We can observe that the energy is smaller everywhere in the GDL system, except 
at intermediate and large heights. Remarkably, the energy of the GDL shock wave is 
lower than the JR after an impact with the wall, however the remnants persist for long 
at larger heights. The MD profile indicates a higher energy at the bottom after an 
impact (c), as compared with both GDL and JR results, but specially with the latter. 
The GDL shock wave is very much damped. It also shows that the impact with the 
bottom wall occurs effectively later, as pointed out when discussing the density profiles. 
In addition, the MD profile shows that the energy vanishes quicker than in the GDL 
solution. Let us examine then the temperature field. 

The most striking difference between the GDL and JR solutions is the temperature 
field, Fig. [6] At large heights, the GDL temperature is one order of magnitude larger 
than the JR. Moreover, the GDL temperature gradient is positive at middle heights (it 
starts to grow) whereas there the JR, like the MD temperature gradient, is negative 
once the shock wave is dissipated. It is clear that the term /xVn helps to sustain large 
temperature gradients in the system, transferring heat from the dense to the dilute 
regions at the top wall. This term is the genuine contribution of the inelastic nature of 
the granular gas to the transport coefficients, although we find no hint in the obtained 
MD profile that the temperature gradient should be positive instead of negative when 
ascending from the dense to the dilute region. As mentioned in the Introduction, the 
presence of the coefficient /i in the heat flux is an exact result of the inelastic Enskog 
equation and the JR approximation fails in describing this new effect. In addition, 
the existence of this term in the heat flux has been already confirmed by computer 
simulation results [25] . 

However one must recall that beyond 40 diameters in height the material gets more 
and more rarefied (Fig. H]) and goes from densities of the order of 1% in packing fraction 
at 40 diameters to about l%o at 60, as obtained from MD results. Therefore one should 
find Knudsen layers when approaching a virtual top wall -in our MD simulations there is 
none, making our hydrodynamic simulations meaningless there. Note on the other hand 
that the temperature field T displayed in Fig. [6j when scaled with the mean free path 
as the relevant unit length, will be proportional to the quantity (f)T displayed in Fig. |5] 
In the latter one can appreciate that the mismatch between JR and GDL is reduced, 
although it still persists. Also, by comparing the three figures (Figs. HI [5] and [6]) that 
the growth of the temperature starts at intermediate heights, when the density is not 
specially small. For this reason, one can conclude that the growth of the temperature 
is a true result of the GDL approach and not a negligible product. 

As the GDL temperature is higher than the JR temperature at the top, the GDL 
solution is more diffusive. Figure |8] shows the vertical component of the heat flux as 
a function of height, where this effect is shown: note the enhanced heat transport at 
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Figure 8: Vertical component of the (reduced) heat flux as a function of height (in units 
of a) at selected times over two oscillation periods, for the JR and GDL simulations. 

For time evolution of the profiles see |38j . 
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intermediate heights, as compared with the JR solution. Unlike the JR case, the GDL 
heat flux consists of two terms, the one coming from the temperature gradient, and 
the one associated, through the coefficient /i, to the density gradient. An analysis of 
the data reveals that both terms have generally opposite signs. The role of the latter 
contribution is to transfer heat from the dense towards the dilute regions at the top, 
while the former brings energy into the granulate, from the high temperature regions at 
the top. Both terms are relevant and contribute in the same order of magnitude. So, 
the heat transfer dynamics is quite different in the GDL and the JR approaches, not 
only at the top but also at the bottom plate when the impacts occur, in such a way that 
gives rise to entirely different solutions for the temperature field. 

In general, the GDL system is less diffusive very close to the plate and more at 
intermediate heights and at the top, as compared with the JR system. The viscosities 
and the cooling term (see Fig. [9]) also follow this pattern. The analysis of the results 
allows us to conclude that in the JR system, most of the energy is dissipated very close 
to the plate, whereas much less is diffused; in the GDL, comparatively, there is less 
dissipation at the plate and more diffusion. 

3.3. Kinetic energy and Mach number 

Figure [10 shows the scaled kinetic energy profiles. An examination of the entire sequence 
shows that the maximum of the kinetic energy is achieved at t = 0.38/" 1 in the GDL 
simulation, at t = 0.42/" 1 in the JR and at t = 0.54/" 1 in the MD. The GDL peak is 
the highest, more than 4 times bigger than the MD, and about 50% bigger than the JR. 
This shows that the GDL solution for the velocity field is also quantitatively different 
from the JR, a consequence of the inelasticity contributing to the viscosities. Leaving 
aside the mismatch at the maximum, the JR and GDL solutions go close to each other, 
and differently from the MD profile, due to the delayed landing of the granular layer in 
the MD simulation. In any case the comparison of the kinetic energy profiles reinforces 
the quite unexpected result that the GDL solution is not closer to the MD, but even 
further away, than the JR. 

Since the GDL temperature is about one order of magnitude higher than JR in 
the dilute region, the Mach number is also smaller. In Fig. CD] we can see how the 
differences are very relevant during the stages (c)-(d), when the layer has achieved its 
maximal extension, and where the JR Mach number is about twice that of the GDL. 
This is another fact showing that the GDL system is more diffusive than the JR. 

The MD curve for the Mach number has been produced using the averaged density 
and temperature fields into Eq. fj29|) for the sound speed, supplied with the equation of 
state (fT3l). 

Unlike JR and GDL approximations, the second MD peak in the Mach number 
is higher than the first one. Anyway GDL predicts better the behavior of the Mach 
number than JR. The values of the Mach number have been computed at the heights 
shown by the red curves in Fig. [TTJ They correspond to the first point, going from 
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Figure 9: The profiles of the cooling term (nT as a function of height (in units of a), at 
selected times over two oscillation periods, for the JR and GDL simulations. Note the 
change in the vertical scale in subfigures (b) and (c). For time evolution of the profiles 

see 138 1 . 
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the dense to the dilute phase, where the packing fraction is 0.1. There we also find 
discrepancies when comparing the MD results with those of JR and GDL simulations. 
This is a consequence of the discrepancies in the density field discussed above. 

4. Conclusions 

In this paper, we have compared the predictions of the Navier-Stokes hydrodynamic 
equations of two dimensional granular gases with MD simulations in a highly nonlinear, 
far-from-equilibrium problem such as the periodic impact of a horizontal piston which 
gives rise to the characteristic pattern formation of the Faraday instability. Given 
that the corresponding Navier-Stokes transport coefficients are not exactly known, two 
different approaches to them have been considered: the JR and GDL approximations. 
While the first approach applies for nearly elastic particles (in fact their forms for the 
transport coefficients are the same as the elastic ones), the latter approximation is much 
more accurate for granular gases (as verified for instance in Ref. [23] by comparing the 
GDL theory with computer simulations at quite extreme values of dissipation) since it 
incorporates the effect of inelasticity on the transport coefficients. In particular, while 
the JR theory neglects the term -/xVn in the heat flux, the new transport coefficient 
fi is clearly different from zero in the GDL theory (see fourth panel of Fig. [T|). After 
comparing both approaches with coarse-grained MD results, we can conclude on the 
following relevant aspects. 

First, we conclude that the granular Navier-Stokes hydrodynamics with the proper 
GDL forms for the transport coefficients 77, 7, k and fi is not capable of reducing 
the discrepancy between discrete particle simulations and hydrodynamic simulations of 
moderately dense, inelastic gases. This quantitative disagreement can be due to the fact 
that while the Navier-Stokes constitutive equations (jl]) and ((5]) for the pressure tensor 
and the heat flux, respectively, apply to first order in the spatial gradients, the problem 
analyzed here might be outside the strict validity of the Navier-Stokes approximation as 
the comparison with MD simulations shows. Surprisingly, the discrepancies between 
theory and simulations decrease if one considers the elastic forms of the transport 
coefficients. We think that there are no physical reasons behind this improvement. 
A similar conclusion has been found in the simple shear flow problem for dilute 
gases since the non-Newtonian shear viscosity rjs(a) to be plugged into the Navier- 
Stokes hydrodynamic equations is better modeled by the elastic shear viscosity than its 
corresponding inelastic version ?7qdl(«) (see Fig. 1 of Ref. [32] where while 7?s decreases 
as decreasing a the opposite happens for t/qdl)- 

Apart from the different ct-dependence of 77, 7 and k, the main difference between 
the JR and GDL approaches lies in the presence of the coefficient fi in the heat flux. 
This new transport coefficient, characteristic of inelastic gases and thus vanishing in 
the JR theory, constitutes the significant contribution to an enhanced heat transfer 
mechanism which leads to a high temperature solution in the dilute region, which is not 
supported by the particle simulations. Even in the dilute region at the top of the system 
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Figure 10: The scaled kinetic energy profiles as a function of height (in units of a), at 
selected times over two oscillation periods for the MD, JR and GDL systems. For time 

evolution of the profiles see [38] . 
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Figure 11: (color online) Mach number for MD system and JR and GDL theories as a 
function of time, along two periods (f^ 1 ) of oscillation of the plate. The red curves 
indicate the variable height (in diameters) which corresponds to the Mach number 
shown. These heights are found as the first point in vertical direction from the plate 

where the packing fraction is 0.1. 



(Fig. [T]), where shock waves should be completely damped, the GDL model produces an 
unrealistic energy excess. There, at densities of the order of 0.001 in packing fraction, 
the coefficient \x is very different from zero. 

Thus, we can conclude that the transport coefficient \x is clearly overestimated 
by the Navier-Stokes approximation and consequently, the influence of the diffusion 
term — /iVn on the heat flux is larger than the one observed in the simulations. As 
mentioned before, these discrepancies do not imply that the GDL approach is deficient 
in any respect. Rather differently, they show the limits of the Navier-Stokes description 
applied to certain regimes of complex granular flow. 

As a matter of fact, in the theory of granular gases it is well accepted that in 
certain cases the Navier-Stokes description is insufficient. It was clearly evidenced on 
the mathematical level, e.g., in jl2], that Burnett-order terms are important for the 
kinetics of granular gases. These terms in the constitutive relations are of second-order 
in the gradients and therefore beyond the Navier-Stokes level of description. In [12] 
it was shown that these terms are even necessary for a consistent description due to 
the lack of scale separation in granular gases. The presence of large gradients is quite 
usual in granular flow, where physical variables may change several orders of magnitude 
within a distance of a few grains, due to inelastic interactions. That typically manifests 
into strong shock waves propagating into the system from the boundaries, where the 
energy source is located. In this way, granular flow is very often supersonic or even 
hypersonic, and in this regime of extreme gradients the Navier-Stokes description may 



25 



reveal inadequate. 

The inclusion of higher-order terms (beyond the Navier-Stokes domain) in the 
constitutive equations for the momentum and heat fluxes might prove a better 
approximation to problems like this one, where the first order in the gradients expansion 
looks insufficient. However, the determination of these nonlinear contributions to the 
fluxes becomes a very hard task if one starts from the revised Enskog equation. In these 
cases it is useful to consider kinetic models with the same qualitative features as the 
Enskog equation but with a mathematically simpler structure [13]. The use of these 
models allows to derive explicit forms for generalized constitutive equations in complex 
states driven far from equilibrium, such as the simple shear flow state jB]. 

In spite of the discrepancies found here, the Navier-Stokes approximation with the 
GDL forms for the transport coefficients is still appropriate and accurate for a wide 
class of flows. Some examples include applications of Navier-Stokes hydrodynamics to 
symmetry breaking and density /temperature profiles in vibrated gases [261 [27], binary 
mixtures jlSj and supersonic flow past a wedge in real experiments pQ, SHI H7] . Another 
group refers to spatial perturbations of the homogeneous cooling state for an isolated 
system where MD results of the critical length for the onsets of vortices and clusters 
[4*8| H9] are successfully compared with the predictions from linear stability analysis [50] 
performed on the basis of the GDL transport coefficients. 

As a summary, the Navier-Stokes theory has shown limitations when exploring the 
highly nonlinear problem of the granular Faraday instability. In particular, the presence 
of rarefied regions where strong transient shock wave fronts propagate seem to justify 
the inclusion of higher order gradients in the transport equations, going beyond the 
Navier-Stokes approximation [42] . In spite of that, both GDL and JR models work 
quite well, although here the main discrepancy is attributed to the term in the heat flux 
coupled to the density gradient, which is the missing contribution in the JR approach 
that the inelastic theory comes to fix. More work has to be done in this respect to 
clarify the conditions under which the Navier-Stokes approximation fails to describe 
appropriately the granular heat transport. Finally, as a complementary route to the 
Navier-Stokes approximation, one could numerically solve the Enskog equation via the 
direct simulation Monte Carlo method [5TJ [52]. Presumably, the numerical solution 
would give better quantitative agreement with MD simulations than the Navier-Stokes 
results reported here. This is quite an interesting problem to be addressed in the 
future for the Faraday and other different problems such as the simple shear flow or 
homogeneous states. 
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